#!/usr/bin/perl -w
use strict;
use FindBin;
use lib ("$FindBin::Bin/..", "/net/cpp-group/Leo/bin");
use chromosomes;
use DBI;
use DBD::Pg;
use linewrap;
use Spreadsheet::WriteExcel::Big;
use xls_formats;

(my $my_dir = $FindBin::Bin) =~ s/(.*)\/$/$1/;

print STDERR <<"HEADLINE";
888888888888888888888888888888888888888888888888888888888888888888888

    xls_protocols_venn_align_kaks

	Description:

888888888888888888888888888888888888888888888888888888888888888888888

HEADLINE

#-----------------------------------------------------------------------------------------
#	Get DB parameters
#
use db_parameters;
use ortho_species;
#my @db_param = get_db_parameters();
#
#	Get DB parameters
#-----------------------------------------------------------------------------------------




#_________________________________________________________________________________________

#	enumerate_intra_clade_pairs

#		prepare a list of pairwise protein comparisons, described by chromosome and 
#			ortholog type

#_________________________________________________________________________________________
sub enumerate_intra_clade_pairs($)
{
	my ($dbh) = @_;
	my $sql_cmd = <<"PL/SQLCMD";
	SELECT 
			gl1.chromosome AS chromosome1, 
			gl2.chromosome AS chromosome2,
			os1.species <> os2.species AS different_species,
			os1.prot_id AS prot_id1, 
			os2.prot_id AS prot_id2, 
			os1.ortholog_type as ortholog_type1, 
			os2.ortholog_type as ortholog_type2 
		INTO 
			TEMP t_intra_clade_pairs
		FROM 
			ortholog_sets os1 JOIN 
			ortholog_sets os2 USING (protocol_id, ortholog_id), 
			taxon.gene_loci gl1, 
			taxon.gene_loci gl2 
		WHERE 
			gl1.gene_id = os1.gene_id AND 
			gl2.gene_id = os2.gene_id AND 
			protocol_id = $curr_protocol_id;
PL/SQLCMD
	my $data = $dbh->do($sql_cmd) or die die $dbh->errstr.".";


	
$sql_cmd =<<"PL/SQLCMD";

	-- orphans other species closest
	INSERT INTO t_intra_clade_pairs
			(chromosome1, chromosome2, different_species, 
			 prot_id1, prot_id2, ortholog_type1, ortholog_type2)
		SELECT 
				gl1.chromosome AS chromosome1, 
				gl2.chromosome AS chromosome2,
				1,
				prot_id AS prot_id1, 
				other_species_prot_id AS prot_id2,
				'orphan ' || species,
				'orphan ' || species,
			FROM 
				predicted_type JOIN 
				closest_other_species cos USING (gene_id, protocol_id) , 
				taxon.gene_loci gl1, 
				taxon.gene_loci gl2 
			WHERE
				gl1.gene_id = cos.gene_id AND 
				gl2.gene_id = cos.other_species_gene_id AND 
				gene_type IN ( 'Orphan singleton inconsistent phylogeny', 
							   'Orphan low dS homolog within species', 
							   'Orphan homolog other species', 
							   'Orphan No homolog', 
							   'Orphan homolog within species', 
							   'Orphan low dS homolog other species' 
							 ) AND 
				protocol_id = $curr_protocol_id;

	-- orphans same species closest
	INSERT INTO t_intra_clade_pairs
			(chromosome1, chromosome2, different_species, 
			 prot_id1, prot_id2, ortholog_type1, ortholog_type2)
		SELECT 
				gl1.chromosome AS chromosome1, 
				gl2.chromosome AS chromosome2,
				0,
				prot_id AS prot_id1, 
				same_species_prot_id AS prot_id2,
				'orphan ' || species,
				'orphan ' || species,
			FROM 
				predicted_type JOIN 
				closest_same_species css USING (gene_id, protocol_id) , 
				taxon.gene_loci gl1, 
				taxon.gene_loci gl2 
			WHERE
				gl1.gene_id = css.gene_id AND 
				gl2.gene_id = css.same_species_gene_id AND 
				gene_type IN ( 'Orphan singleton inconsistent phylogeny', 
							   'Orphan low dS homolog within species', 
							   'Orphan homolog other species', 
							   'Orphan No homolog', 
							   'Orphan homolog within species', 
							   'Orphan low dS homolog other species' 
							 ) AND 
				protocol_id = $curr_protocol_id;

	CREATE INDEX i__t_intra_clade_pairs ON t_intra_clade_pairs (ortholog_type1, ortholog_type2, different_species, prot_id1, prot_id2);
	CREATE INDEX i__t_intra_clade_pairs1 ON t_intra_clade_pairs (ortholog_type1, ortholog_type2, different_species, prot_id1, prot_id2) WHERE ortholog_type1 ~ 'to' AND ortholog_type2 ~ 'to';
PL/SQLCMD
	$dbh->do($sql_cmd) or die DBI::errstr."\n$sql_cmd.";
	
}




#_________________________________________________________________________________________

#	get_histogram

	# calculate counts, quartiles and percentile histograms
	# output format is

	#			counts
	#			quartiles
	#			mean
	#			percentiles

#_________________________________________________________________________________________
sub get_histogram(@)
{
	my @data = sort {$a <=> $b} grep {defined $_} @_;

	# calculate mean
	my $sum = 0.0; 
	$sum += $_ for (@data);
	my $mean = 0.0;
	$mean = $sum / @data if @data;
	
	my $count				= scalar @data;
	my $count_undefined		= scalar grep {!defined $_} @_;
	my $count_defined		= scalar @data;
	my @indices = (int($count_defined * 0.25), int($count_defined * 0.50), int($count_defined * 0.75));
	my @quartiles = @data[@indices];
	
	@indices = (0..100);
	$_ = int($_ * $count_defined / 100) for @indices;
	$indices[-1] = @data - 1;
	$indices[0] = 0;
	my @percentiles = @data[@indices];

	return [[$count, $count_defined, $count_undefined], [@quartiles], $mean, [@percentiles]];
}













#_________________________________________________________________________________________

#	do_select_subset

#		get count or histogram of data divided by all, protocol1 only , protocol2 only, 
#				common or all protocol1, all protocol2

#_________________________________________________________________________________________
sub do_select_subset($$$$$)
{
	# new
	my ($dbh, $what_field, $sql_filter, $count, $return_data) = @_;



	my $left = $count ? " LEFT " : "";
	my $total_id_filter = $id1 + $id2;
	my $sql_cmd = <<"PL/SQLCMD";
    SELECT
            $what_field
        FROM 
			t_intra_clade_pairs LEFT JOIN
            orthologs.align_kaks aks USING (prot_id1, prot_id2)
PL/SQLCMD
	$return_data = $dbh->selectall_arrayref($sql_cmd) or die DBI::errstr."\n$sql_cmd.";
	$sql_cmd =~ s/[ \n]+/ /g;
	#print STDERR "\t", $sql_cmd, "\n";
	#print STDERR "Returned ", scalar @$data, " from panda\n";
}


#_________________________________________________________________________________________

#	get_ortholog_histograms

#		get count or histogram of data for orthologs, 1:1 orthologues, manys, paralogues

#_________________________________________________________________________________________

pseudos
pseudos against non-pseudos

orphans
closest other species
closest same species


sub get_ortholog_histograms($$\@\@)
{
	
	my ($dbh, 
		$what_field,
		$return_data, # return data
		$column_names,
        ) = @_;
	print STDERR "get_ortholog_histograms\n";

	#
	#	all orthologs
	#
	print STDERR "\tGet all orthologues data\n";
	{
		my $data;
		do_select_subset($dbh, 
						 $what_field,
						 "WHERE ortholog_type1 ~ 'to' AND ortholog_type2 ~ 'to' AND different_species",
						 $data);
		my @data = grep {$_->[0];} @$data;
		push(@$return_data, get_histogram(@data));
		push(@$column_names, "All orthologues");
	}




	#
	#	pseudos
	#
	# select manys orthologues
	{
		my $data;
		do_select_subset($dbh, 
						 $what_field,
						 "WHERE (ortholog_type1 ~ 'to' OR ortholog_type1 ~ 'inparalog') ".
						 "AND ortholog_type1 = ortholog_type2 AND same_species",
						 $data);
		my @data = grep {$_->[0];} @$data;
		push(@$return_data, get_histogram(@data));
		push(@$column_names, "All inparalogues");
	}
	

	#
	#	per ortholog type
	#
	my %data;
	for my $ortholog_type("1 to 1", "1 to m", "m to 1", "m to m")
	{
		print STDERR "\tGet $ortholog_type orthologues data\n";
		do_select_subset($dbh, 
						 $what_field,
						 "WHERE ortholog_type1 = '$ortholog_type' and ".
						 "ortholog_type2 = '$ortholog_type'  AND different_species",
						 $data{$ortholog_type});
		my @data = grep {$_->[0];} @{$data{$ortholog_type}};
		push(@$return_data, get_histogram(@data));
		push(@$column_names, "$ortholog_type orthologues");
	}


	# select manys orthologues
	{
		print STDERR "\tGet $ortho_name1 manys orthologues data\n";
		my @manys = (@{$data{'m to 1'}}, @{$data{'m to m'}});
		push(@$return_data, get_histogram(@manys));						
		push(@$column_names, "$ortho_name1 manys orthologues");

		print STDERR "\tGet $ortho_name2 manys orthologues data\n";
		my @manys = (@{$data{'1 to m'}}, @{$data{'m to m'}});

		my @data = grep {$_->[0];} @manys;
		push(@$return_data, get_histogram(@data));
		push(@$column_names, "$ortho_name2 manys orthologues");
	}


	# select paralogues
	print STDERR "\tGet all in-paralogue data\n";
	{
		my $data;
		do_select_subset($dbh, 
						 $what_field,
						 "WHERE (ortholog_type1 ~ 'to' OR ortholog_type1 ~ 'inparalog') ".
						 "AND ortholog_type1 = ortholog_type2 AND same_species",
						 $data);
		my @data = grep {$_->[0];} @$data;
		push(@$return_data, get_histogram(@data));
		push(@$column_names, "All inparalogues");
	}

	# select Non-orphan paralogues
	print STDERR "\tGet Non-orphaned in-paralogue data\n";
	{
		my $data;
		do_select_subset($dbh, 
						 $what_field,
						 "WHERE ortholog_type1 ~ 'to' ".
						 "AND ortholog_type1 = ortholog_type2 AND same_species",
						 $data);
		my @data = grep {$_->[0];} @$data;
		push(@$return_data, get_histogram(@data));
		push(@$column_names, "Non-orphaned inparalogues");
	}

	# select Orphan paralogues
	print STDERR "\tGet all orphaned in-paralogue data\n";
	{
		my $data;
		do_select_subset($dbh, 
						 $what_field,
						 "WHERE ortholog_type1 ~ 'inparalog' ".
						 "AND ortholog_type1 = ortholog_type2 AND same_species",
						 $data);
		my @data = grep {$_->[0];} @$data;
		push(@$return_data, get_histogram(@data));
		push(@$column_names, "Orphaned inparalogues");
	}
}

#_________________________________________________________________________________________

#	write_range_data_to_xls

#		send histogram etc to xls

#_________________________________________________________________________________________
sub write_range_data_to_xls($$$\@\@\@$)
{
	my ($worksheet, $format, $cnt_blocks, $block_names, $row_names, $data, $data_format) = @_;

	# write data
	{
		my $col = 1;
		for $data(@$data)
		{
			$worksheet -> write_col(2, $col, $data->[0]);
			$worksheet -> write_col(4, $col, $data->[1], $data_format);
			++$col;
		}
	}
	
	# write percentages (Y column) every 8
	my @range = map {$_ / 100} (0..100);
	$worksheet->write_col(10, $_ * 8 - 1, \@range, $format->{blue_bold_percent}) 
																for (1.. $cnt_blocks);

	# write block names
	$worksheet->write(0, $_ * 8 + 1, $block_names->[$_], $format->{red_bold_yellow}) 
															for (0.. ($cnt_blocks - 1));
	# write header row and column
	{
		my $col = 0;
		for (@$row_names)
		{
			if ($_)
			{
				$worksheet->write(1, $col, $_, $format->{blue_bold});
				$worksheet->write(9, $col, $_, $format->{blue_bold});
			}
			++$col;
		}
	}
	$worksheet -> write_col(2, 0, 
							['count', 'count missing', '="25%"', "median", '="75%"', 'mean'], 
							$format->{blue_bold});
}


#_________________________________________________________________________________________

#	get_specific_histogram_for_orthologs

#		send histogram etc to xls

#_________________________________________________________________________________________
sub get_specific_histogram_for_orthologs($$$$\%$)
{
	my ($dbh, $workbook, $what_field, $name, $format, $format_for_data) = @_;

	print STDERR "Get $name histogram...\n";
	my (@results, @column_names);
	get_ortholog_histograms($dbh, $what_field, @results, @names);

	my $worksheet = $workbook->addworksheet($name);
	write_range_data_to_xls($worksheet, $format, $format_for_data, 
							@column_names, @results);

}





#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
#
#	Get data from panda
#
#
my $dbh = connect_to_panda();
$dbh->{RaiseError} = 0;

system ("mkdir -p $dir_orthologs_output/excel") and die;
my $workbook  = Spreadsheet::WriteExcel::Big->new("$dir_orthologs_output/excel/orthologs.align_kaks.xls");
die "Problems creating new Excel file: $!" unless defined $workbook;

my %formats;
get_xls_formats($workbook, %formats);

$dbh->do ("set work_mem=160384") or die DBI::errstr.".";


	#
	# get protocol_ids and comparison pairs
	#
	my $curr_protocol_name;

	get_specific_histogram_for_orthologs($dbh, $wookbook, 'dnds', 'dN/dS', %formats, 'double_digit');
	get_specific_histogram_for_orthologs($dbh, $wookbook, 'ds', 'dS', %formats, 'double_digit');
	get_specific_histogram_for_orthologs($dbh, $wookbook, 'dn', 'dN', %formats, 'double_digit');
	get_specific_histogram_for_orthologs($dbh, $wookbook, 'percent_id', '%id', %formats, 'percent');
	get_specific_histogram_for_orthologs($dbh, $wookbook, 'align_len', 'aligned residues', %formats, undef);
	get_specific_histogram_for_orthologs($dbh, $wookbook, 'kappa', 'kappa', %formats, 'double_digit');
	get_specific_histogram_for_orthologs($dbh, $wookbook, 'tau', 'tau', %formats, 'double_digit');
	get_specific_histogram_for_orthologs($dbh, $wookbook, 'coverage', 'coverage', %formats, 'percent');
	get_specific_histogram_for_orthologs($dbh, $wookbook, 'aligned_residues / less(len1, len2)', 'shorter_coverage', %formats, 'percent');
	get_specific_histogram_for_orthologs($dbh, $wookbook, 'aligned_residues / more(len1, len2)', 'longer_coverage', %formats, 'percent');
	get_specific_histogram_for_orthologs($dbh, $wookbook, 'count_s / aligned_residues', '% synonmous sites', %formats, 'percent');
	get_specific_histogram_for_orthologs($dbh, $wookbook, 'more(len1, len2) - less(len1,len2)', '% synonmous sites', %formats, undef);
	

$workbook->close();
print STDERR "Finished.\n";


Orphans to homolog

genomic length
exon numbers
coding length




PL/SQLCMD
	$dbh->do($sql_cmd) or die DBI::errstr."\n$sql_cmd.";
	my $sql_cmd = <<"PL/SQLCMD";
PL/SQLCMD
	$dbh->do($sql_cmd) or die DBI::errstr."\n$sql_cmd.";



